A Prototype of Dynamical Quantum Phase Transitions of Glassy Systems 
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We explore a typical behavior of glassy systems driven by quantum dynamics in a quenching pro- 
tocol. Using the random energy model in a transverse field, we calculate several types of dynamical 
quantum amplitude. The result shows a freezing transition at some critical time, indicating the 
dynamical quantum phase transition. The behavior is understood by the partition-function zeros 
in the complex temperature plane. We discuss the properties of the freezing phase as a dynamical 
chaotic phase, which are contrasted to those of the spin-glass phase in the static system. 



Phase transitions are one of the most fascinating top- 
ics in statistical physics. After the long investigations, 
equilibrium phase transitions are clarified as singulari- 
ties of free energies and become known to be treatable 
in the current framework of statistical mechanics [H-dj]. 
To analyze those singularities, sophisticated techniques 
and concepts, such as the renormalization group theory 
and the scaling hypothesis, have been founded and nowa- 
days equilibrium phase transitions of classical systems are 
fairly well understood. On the other hand, phase transi- 
tions or dynamical singularities of quantum systems are 
far less understood, though rising interests in quantum 
computation require better comprehension of those phe- 
nomena. Recent experimental developments enable to 
observe time evolutions of quantum many-body systems 
in detail which also demands better theoretical 

clarifications of dynamical behavior of quantum systems. 

Recently, dynamical singularities of quantum systems 
in a typical experimental set up, the so-called quantum 
quench, are focused 0-0- In the set up, a quantum sys- 
tem, whose Hamiltonian with a parameter Y is written 
as "H(r), is prepared in a quantum state |V>o) dominated 
by ri(Yo), and the parameter is suddenly changed as 
ro — > Yi at t — 0. After the quench, the wave func- 
tion of the system at time t > is given by 

-* tW(ri U). 
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Due to the sudden quench, the behavior of this wave func- 
tion can be drastically different from the one in adiabatic 
process. It is actually reported that another energy scale 
different from the adiabatic one exists and dominates the 
system, which can lead to novel non-analyticity in dy- 
namical behavior [7- 9]. Interestingly, those dynamical 
singularities can be analyzed by the theory of zeros of 
equilibrium partition function [H El [t| . This provides a 
possibility to analyze a wide range of dynamical singu- 
larities of quantum systems, with the current maturity 
of the theory of zeros. 

Based on this possibility, in this paper, we provide a 
prototype of dynamical singularities of glassy systems. 
For this, we focus on the so-called random energy model 
(REM) [HOII in a transverse field (TREM) introducing 
quantum nature [l2|-[T5| . Despite its simplicity, the REM 



is known to show many nontrivial properties of glassy 
systems and is treated as an ideal platform to study sys- 
tems with many metastable states [l6l . [l7| . which nat- 
urally motivates us to firstly investigate the TREM to 
examine dynamical singularities of glassy systems. In 
addition, the distribution of the partition- function zeros 
is well known in the REM [l5|, [l8| , and we can expect 
that the present analysis reveals a clear relation between 
the dynamical singularities and zeros. 

The Hamiltonian of the TREM is defined as the limit of 
the fully-connected Ising spin-glasses (SGs) with p-body 
interactions rl — H v 



H g where 



N 
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a x,z are p am j matrices, N is the number of spins, and the 
interaction J^.-.i is independent identically distributed 
(i.i.d.) from Gaussian with zero mean and variance 
p\/2N p - 1 . The REM is produced by taking the limit 
linip^oo H p = #rem EH- Energy levels of if rem lose 
the correlations among each other in the limit p — > oo, 
and thus energy levels of the REM Hamiltonian -//rem 
become i.i.d. Gaussian variables with zero mean and 
variance N/2. The TREM Hamiltonian is given by 
%trem = -Hrem + H g , and hence in the a z represen- 
tation the diagonal elements are random Gaussian vari- 
ables coming from -Srem, and the off-diagonal ones arise 
from H g . 

The TREM has three thermodynamic phases [l2| ■ One 
is the classical paramagnetic (CP) phase which dom- 
inates the system for small Y and large temperature 
T = and another one is the SG phase appearing at 
low T and small Y. The last one is the quantum param- 
agnetic (QP) phase emerging for large Y. The free energy 
densities are given by 



In 2 



for CP 
-Vhrir for SG 

-iln(2cosh/3r) for QP 
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First-order phase transitions from QP to CP or SG 
phases occur at critical values T C (T). Clearly, the CP and 
SG phases do not depend on T, and also the QP phase is 
not influenced by -Srem since / is completely identical to 
— (1/./V/3) InTr exp(— (3H g ). These imply the strong sim- 
plicity of the TREM. Actually, Jorg et al. showed that 
perturbations from T = and T = oo do not change the 
eigenfunctions in both two cases in the thermodynamic 
limit [14j , which can be schematically written as 



Tr e 



Tle -PH g f or r>r c (T) 
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This property makes the following analysis extremely 
simple. 

Hereafter, we consider a quantum quench protocol 
from large to small T, and examine the dynamical 
singularities in the TREM. Following Ref. [9|, we in- 
vestigate the return amplitude G(t) = (ipo\tp(t)) = 
(-00 1 exp(~ itH)\ipo), especially focusing on the ground 
state for large T. Namely, we put |^o) = IQP) where 
|QP) denotes the ground state of H g which includes all 
the cr z -basis with an equal weight 1 / V2 N . 

Thanks to the specialty of the TREM ([5]) , we can eas- 
ily say |G(i)| 2 = 1 if the quench does not go across the 
critical value r c (0). On the other hand, when the quench 
goes across r c , |G(i)| 2 is expected to exponentially de- 
crease as time t grows, i.e. |G(i)| 2 oc exp(— Ng(tj). To 
assess the rate function g(t), Eq. ([5]) again makes the 
problem quite simple. Using this, we can rewrite 



G(t) = (QP| ( 



|QP) = -Lz RE m(^), (6) 



where we write the partition function of the REM under 
the inverse temperature j5 as Zrem(/3). Hence, g(t) is 
directly related to the free energy of the REM under the 
imaginary temperature. 

The calculation of Eq. (J6j is a simple task in 
the REM. Using the spectral representation G(t) — 

2 Sn=i ex P( — iEnt), we divide the square of the am- 
plitude |G(i)| 2 into time-dependent and -independent 
parts. Taking the random average, we obtain 



[|G(t)| 2 ] - exp(-AHn2) + exp(-A^ 2 /2). 



(7) 



We naively expect that the first term is much smaller 
than the second one since the former involves 2 JV -terms 
and the latter 2^(2^ - 1) - 2 2N . However, when 
t becomes large enough, the second Gaussian term is 
considerably smaller than the first term and the time- 
independent part becomes important. Thus, we have a 
dynamical phase transition. The rate function is given 
by 



^ for t < t c = V21n2 
In 2 for t > t c 



(8) 



We see that the amplitude freezes at t 



To obtain further information about \ip(t)), overlaps 
with other wave functions are also useful. As such ref- 
erence wave functions, we here introduce wave functions 
reflecting finite temperature properties 



-/3 //rem/2 



IQP) 



V(QP|e-^™|QP) 
and consider the overlap with \ip(t)) 



g 1/2 (/m) = (f3\m) = 



^REM 



l+it 



(9) 
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We can easily confirm that the quantum-mechanical av- 
erage of physical quantities commutative with -Hrem by 
\(3) completely agrees with the statistical mechanical av- 
erage at inverse temperature (3 [TT^, [2(| . Concerning non- 
commutative quantities, correct expectation values can 
also be derived by appropriately extending the Hilbert 
space of |/3) [2l| - |23j or constructing a certain ensem- 
ble [H, |25j. These considerations support the use of \j3) 
for capturing properties at finite temperatures. 

To treat G(t) and Gi/ 2 (/3,t) in the same manner, we 
introduce a generalized overlap with index k as 



^REM 



(k/3 + it) 



(2")i-* (ZremO?))* 



(11) 



This reduces to Eq. (HUJ) at k = 1/2, and G k (0,t) = 
G (l3,t) = G{t). Furthermore, Gi(/3,t) denotes 
the thermal average of the time evolution operator 
exp(— HHjiem) and is worth studying. 

To handle the average [|Gfc| 2 ], we use the replica 
method 



1 



lim 



(2^)2(1-^') n "^i n2 ™ 2k 
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REM 
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We treat n\ and n 2 as integers and consider the ana- 
lytic continuation to real variables after the calculation. 
The calculation goes along the standard manner [15| and 
the result is expressed in terms of the order parame- 
ter matrix Q with the size ri\ + n\ + n 2 . We obtain 
[|G fe | 2 ] = lim„ 1 ^ilim„ 2 ^_ 2 fcCxp(-7Vg) w here 



9 = 3 ) Alab 



a , b 
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2 Z.^ V ab 

a,b 



(+-) 

b 



/3(kf3 + it) ( +2 ) P(kp -it) ( _ 2 ) 

Z^ q ab 9 2-1 2-1 q ab 



a b 

-2(1 -fc)ln2-s. 



a b 



(13) 



q = or 1 represents a matrix element of Q and is de- 
termined by the extremized condition. The superscript 
denotes the block in the matrix Q: + (—) is for the first 
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(second) block with the size n%, and 2 is for the third 
with 7i 2 . The subscript denotes replica numbers in each 
block. The number of configurations for a given order 
parameter Q is denoted by exp(TVs) and s plays the role 
of the entropy density. 

Following the discussions in [l5[ , we can find four pos- 
sible saddle-point solutions. In the replica symmetric 
(RS) level, we have two solutions and the rate function 



9k = 
PI: 



lim ni _n lim r , 



-2k 9 is calculated as follows: 



This solution is denoted by the identity matrix Q = 
1. That is: q^' = ' = q$ = 8 a ,b and zero for 
the other elements. The corresponding entropy is 
s = (2ni + 7i 2 ) In 2. Then, the rate function is 



f 



9k 



k(l - k) 



ft 2 



(14) 



P2: The first and second blocks has an identical config- 
uration as q£> = q { ~ ] = qF} = q { ^" ] = 8 a , b [Hj]. 
The entropy is given by s — (rix + ?i 2 ) In 2 and we 
obtain 



(15) 



Other two solutions are obtained by assuming the one- 
step replica symmetry breaking (1RSB). In the present 
problem, we divide 71 1 replicas into n\/m\ blocks of size 
mi, and q a b takes 1 if a and b are in the same block 
and otherwise. The 7i 2 replicas are similarly divided 
into 712/7712 blocks of size 7712- Two different solutions 
emerges according to the value of q^ 

1RSB1: tjf( +_ ) = 0. The entropy is s = (2n 1 /m 1 + 
712/7712) In 2 and the rate function becomes 



gk{mi,m 2 ) = 2(1 - k) In 2 - 



t/ 



21n2 



-7711 



771l 



k/3 2 



-7712 



2fcln2 



7712 



(16) 



Extremization with respect to 7112 give m\ = 
2Vln 2//3 = /3 C / p. On the other hand, the extrem- 
ization of mi gives a physically-unacceptable con- 
dition m\ = p c /y/k 2 P 2 — t 2 . Hence, we get one 
possible solution in this ansatz 

g k (l, 771*) = U 2 - l(kp - P c ) 2 + 2(1 - k) ln2. (17) 

1RSB2: We assume that q^ > takes the same form as 
and q^\ Then, 



^(7711,777-2) = 2(1 -fc)ln2- k I p l m l 

kp 2 2fcln2 
H — jq-m 2 "I ■ 

2 7712 



In 2 

771 1 



(18) 



The extremization gives = p c /2kp and — 
Pc/ P, both of which are acceptable in contrast to 
the 1RSB1 case. Hence, we have three possible 
solutions in this case 

g k (l,m* 2 ) = - Up - +2(l-fc)ln2, (19) 

<7 fe K,l) - ^(/3-/3 c ) 2 + 2(l-£:)ln2, (20) 
g k (m* 1 ,m* 2 ) = 2(1 -k)ln2. (21) 

Comparing the above solutions and considering the 
RSB constraint < mi, 7712 < 1, we specify the dominat- 
ing solutions, to find that the set of solutions depend on 
k. For 1 < k, the solutions are PI, P2, and 1RSB2 with 
(7711,7712) = {m\, 1), (ml,m* 2 ). For 1/2 < k < 1, 1RSB1 
is included to the above four solutions. For k < 1/2, 
1RSB2 with (1, m|) appears instead of that with (m*, 1). 
We show the phase diagram at k = 1 and 1/2 in Fig. [T] 
In PI and 1RSB1 phases, the rate function depends on t 
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FIG. 1. Phase diagram for k — 1 and k — 1/2 on /3-t plane. 



and the amplitude decays in t. Then, in P2 and 1RSB2 
phases, the amplitude ceases to decrease and is frozen to 
a fixed value. 
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FIG. 2. The rate function g(t) in the discrete REM with 
N = 80 and a = 1.3. Each line represents the result from a 
single sample. 




FIG. 3. gi/ 2 (/3,t) in the same calculation as Fig. [2] 



Our result of [|Gfc| 2 ] is compared with the distribution 
of the partition-function zeros in the complex-/? plane. 
The zeros of Z-hem(0 + H) was calculated in Refs. [i&ITsT] 
and a similar phase diagram including the PI, P2, and 
1RSB phases was found. In the present problem, the rate 
function for a given is frozen to a fixed value in the 
P2 and 1RSB2 phases and it seems to be hard to recog- 
nize the difference between these two phases. However, it 
is known that the density of zeros two-dimensionally dis- 
tributes in the P2 phase but no zeros exist in the 1RSB2 
phase [TH, [IH , which implies a qualitative difference. 

Concerning the 1RSB2 phase, it is easy to interpret the 
plateau of gk with respect to t, since this corresponds to 
the SG phase which is well understood. The SG transi- 
tion of the REM occurs as the freezing of the system into 
the ground state. This directly leads to gt in the 1RSB2 



phase since the ground state of the REM, |SG), gives the 
overlap \(SG\^{t))\ 2 = 1/2* 

On the other hand, such a freezing does not occur in 
the P2 phase. Instead, it is expected that the dominant 
states for the free energy rapidly switch as t changes, 
which implies a kind of chaos emerges in the overlap in 
the P2 phase. This speculation is based on the recently- 
proposed relation between the distribution of zeros and 
the chaos effect of SGs (26[. The chaos effect of SGs 
means that the spin configuration drastically changes as 
the physical parameters, such as temperature or external 
field, slightly vary [27| ■ In the mean- field level, the chaos 
effect of SGs can be interpreted as the changes of domi- 
nant thermodynamic pure states [28[ , which produces the 
two-dimensionally distributing zeros [26]. These consid- 
erations naturally lead to the connection between a kind 
of chaos in dynamic behavior and the P2 phase. 

To confirm this assertion, we numerically calculate the 
rate function. From a technical reason, we treat the dis- 
crete version of the REM, where 2 Ar -energy levels are 
denoted by M + 1-discrete binomially-distributed vari- 
ables [1^, |3(|. This model gives a phase transition for 
a = M/N > 1, whose physical nature is the same as the 
REM. We plot the rate function g(t), which is equivalent 
to g\/i{fi — 0,i), at a = 1.3 in Fig. [2j The plotted curves 
are taken from 10 samples of random distributions. We 
can find that for small t, where the PI phase dominates, 
all the samples give an identical smooth curve, while for 
large t singular fluctuations depending on samples appear 
around the plateau of the P2 phase git) = In 2. This sin- 
gular behavior is just the chaotic state mentioned above. 
To make further comparison, in Fig. [3] we also plot the 
data at /3 = 8 where the 1RSB2 phase dominates. Al- 
though the sample dependence is large both for the P2 
and the 1RSB2 cases, there is no singular fluctuations 
for the latter case. This difference originates from the 
partition-function zeros, which are present and absent in 
the P2 and 1RSB2 phases, respectively. 

As a concluding remark, we should notice that this 
chaotic behavior in the P2 phase is not specific to 
the REM. Actually, the existence of the P2 phase is 
confirmed in other SG models (2 6l l3lj such as the 
Sherrington-Kirkpatrick model [32J and the spherical 
model 33]. This fact implies a universal aspect of dy- 
namical singularities of glassy systems, which motivates 
to further examine the dynamical singularities of glassy 
systems. 
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